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We present a systematic study of quantum phases in a one-dimensional spin-polarized Fermi gas. 
Three comparative theoretical methods are used to explore the phase diagram at zero temperature: 
the mean-field theory with either an order parameter in a single-plane- wave form or a self-consistently 
determined order parameter using the Bogoliubov-de Gennes equations, as well as the exact Bethe 
ansatz method. We find that a spatially inhomogeneous Fulde-Ferrell-Larkin-Ovchinnikov phase, 
which lies between the fully paired BCS state and the fully polarized normal state, dominates most 
of the phase diagram of a uniform gas. The phase transition from the BCS state to the Fulde- 
Ferrell-Larkin-Ovchinnikov phase is of second order, and therefore there are no phase separation 
states in one-dimensional homogeneous polarized gases. This is in sharp contrast to the three- 
dimensional situation, where a phase separation regime is predicted to occupy a very large space 
in the phase diagram. We conjecture that the prediction of the dominance of the phase separation 
phases in three dimension could be an artifact of the non-self-consistent mean-field approximation, 
which is heavily used in the study of three-dimensional polarized Fermi gases. We consider also the 
effect of a harmonic trapping potential on the phase diagram, and find that in this case the trap 
generally leads to phase separation, in accord with the experimental observations for a trapped gas 
in three dimension. We finally investigate the local fermionic density of states of the Fulde-Ferrell- 
Larkin-Ovchinnikov ansatz. A two-energy-gap structure is shown up, which could be used as an 
experimental probe of the Fulde-Ferrell-Larkin-Ovchinnikov states. 
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I. INTRODUCTION 

Since the successful demonstration of a magnetic Fes- 
hbach resonance^ and the creation of optical lattices^, 
ultracold atomic Fermi gases have become a topic of great 
current interest^. Thanks to these key tools, the inter- 
atomic interactions and even the dimensionality of ul- 
tracold atomic Fermi gases can be easily tuned, which 
makes them ideal candidates to simulate novel quan- 
tum many-particle systems. Therefore, an intriguing 
opportunity is opened for studying some long-standing 
problems, such as the crossover from Bardeen-Cooper- 
Schrieffer (BCS) superfluidity to Bose-Einstein conden- 
sate (BEC) 0, S, fi (3, H, 3, and models of high tem- 
perature superconductivity. These remarkable prospects 
have attracted attention from many researchers, rang- 
ing from condensed matter physics to atomic molecular 
and optical physics, and even particle and astro physics. 
Experimentally, superfluidity of an ultracold Fermi gas 
at the strongly interacting B CS-BEC crossover has been 
strikingly demonstratedflQ, llil, [13, H Q E E [13, H 
fl9l.[20l.l2l]|. This is a landmark achievement in the history 
of physics. 

Recent experiments have now generated ultracold 
atomic Fermi gases with flnite spin polarizationjll, [H, 
[13, [H, [H, [13 ■ That is, the two spin components have 
unequal populations. However, the physical understand- 
ing of the ground state of a polarized atomic gas remains 
an open question. The standard BCS model - though 
not quantitatively accurate for strong interactions - is 



still qualitatively correct when there is no spin polar- 
ization. This simply involves Cooper pairing between 
spin up and spin down atoms with opposite momenta 
at the same Fermi surface. A polarized Fermi gas can- 
not be explained within standard BCS theory because 
the Fermi surfaces of the two spin components are mis- 
matched. Non-standard forms of pairing must exist to 
support superfluidity in this polarized environment. 

The study of polarized Fermi gases can be traced back 
to the middle of the twentieth century, soon after the 
seminal BCS theory paper. Similar theoretical propos- 
als were independently given by Fulde and Ferrell [2a|, 
and Larkin and Ovchinnikov [29] (FFLO). These authors 
suggested that Cooper pairs may acquire a flnite center- 
of-mass momentum jsO]. In such an ansatz, the two mis- 
matched Fermi surfaces can overlap, thereby supporting 
a spatially inhomogeneous superfluidity. The search for 
the existence of the predicted FFLO state has lasted for 
more than four decades. Only very recently has there 
been indirect experimental evidence for observing such 
states in the heavy fermion superconductor CeCoIns [3l|. 
Due to the shrinkage of the available phase space for pair- 
ing, the FFLO state is now thought to be very fragile in 
three dimensions. Alternative pairing scenarios include: 
Sarma superfluidity [32, .33,,^], a deformed Fermi surface 
[H, [IB, [131 , and breached pairing jsH]. However, at zero 
temperature these phases may suffer from an instability 
towards phase separation. As a result, a phase separa- 
tion regime consisting of a conventional BCS superfluid 
and a normal fluid may be favored in three dimensions 
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The above theoretical issues were not completely re- 
solved in current measurements on polarized ^Li gases 
near a broad Feshbach resonance, carried out at MIT 
[13, [13, [3, m and Rice university [H, 111- Though 
a clear quantum phase transition from a superfluid to 
normal state was observed ^2^, the nature and the 
order of the transition could not be determined due 
to the finite experimental resolution. The presence 
of a harmonic trap in these experiments caused addi- 
tional difficulties in interpreting the experimental re- 
sults. A number of theoretical papers have sought to ex- 
plain these experiments on polarized atomic Fermi gases 
|40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,M, 
56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69,i3l,El|, 
II, 13 , 14,, J^, M, Jl, J&., 19,, M, SI, M, M,M]- > From 
these analyses, the issues that require timely clarification 
may be summarized as follows: 

(A) Structure and detection of FFLO states. De- 
spite a long history, the precise structure of the FFLO 
states remains elusive |30|- Current investigations of 
FFLO states rely mostly on the use of a single-plane- 
wave form for the pairing order parameter A(x), where 
A(x) = Ao exp[iq • x], as initially proposed by Fulde and 
Ferrell (PF). Here q is the center-of-mass momen- 
tum of the Cooper pairs, and the ansatz implies that the 
magnitude of the order parameter and density is con- 
stant in space [3, lil, [HH] ■ The resulting window for the 
FFLO state in parameter space turns out to be very nar- 
row [31 • Can we expect a larger parameter range after 
an optimization of the FFLO proposal? Indeed, by im- 
proving the form of the order parameter to the Larkin 
and Ovchinnikov (LO) type, A(x) c>c cos[q-x], Yoshida 
and Yip have found recently that the FFLO state became 
more stable [H^]. On the other hand, so far there is no 
conclusive evidence for the experimental observations of 
FFLO states 

(B) Intrinsic reason for phase separation. The narrow 
window of the FFLO state may require phase separa- 
tion to fill the gap between BCS and FFLO phases in 
the phase diagram [3§|. Experimentally, a shell struc- 
ture in the density profile of polarized Fermi gases was 
observed [23 . [2^ . suggesting an interior core of a BCS 
superfiuid state with an outer shell of the normal com- 
ponent. Phase separation in trapped systems, however, 
cannot be used as a definitive support of the existence 
of phase separation in a homogeneous gas, since the trap 
favors separation. 

(C) Quantitative approach for polarized Fermi gases 
at the BCS-BEC crossover. A more serious problem is 
the validity of the mean-field approach. The experiments 
were done in the strongly interacting BCS-BEC crossover 
regime, where for the quantitative purpose strong pair 
fiuctuations must be taken into account 0, 0, 0, 13| • -Be- 
cause of the lack of reliable knowledge of the superfiuid 
phase, these pair fiuctuations are usually onl y co nsidered 
above the superfiuid transition temperature |49l . . For 
the same reason, numerical quantum Monte Carlo simu- 



lations have been restricted to the normal state [H^ . [63 | 
and hence cannot provide useful information for the su- 
perfiuid state. 

To gain a qualitative insight into these crucial points, 
in a recent Letter [sH], we have considered a polarized 
Fermi gas in one dimension (ID) at zero temperature. In 
this case the model in free space is exactly soluble via 
a Bethe ansatz solution [H, \m, S [sl, [111, 19^]. We 
have estabhshed the ID phase diagram of the polarized 
gas, both in the uniform situation and in the experimen- 
tally important trapped environment. Complemented by 
a mean-field Bogoliubov-de Gennes (BdG) calculation, 
we have shown that a phase similar to the FFLO-type 
polarized superfiuid is the most widespread in the phase 
diagram. Using a local density approximation to account 
for the harmonic trapping potential, we have found that 
the trap generally leads to phase separation, with at least 
one FFLO-type phase present at the trap center. 

In this paper, we discuss these results in greater detail, 
and compare them to other approximations. We partic- 
ularly focus on the self-consistent BdG method, which 
we previously treated briefiy|i85|]. To address the issue 
of the different possible FFLO structures, we present a 
simplified mean-field calculation with a single-plane- wave 
assumption for the order parameter, and compare it with 
the self-consistent BdG results. These systematic inves- 
tigations give rise to a comprehensive quantitative un- 
derstanding of the ID polarized Fermi gas. We note that 
a qualitative picture was also obtained in earlier some 
works, which were based on a non-perturbative bosoniza- 
tion analysis jo^l or a mean-field approximation with an 
additional assumption on the single-particle energy spec- 
trum (93. [95|. However, the resulting phase diagram was 
not conclusive, and the nature of the transition from BCS 
to FFLO states was under debate [i^l- 

Strictly speaking, any mean-field approach is only valid 
in the weak coupling limit. As the interaction strength 
increases, the pair fiuctuations become increasingly im- 
portant, and therefore have to be taken into account. 
This is particularly noticeable in ID, where true long- 
range order is completely destroyed by fiuctuations in a 
homogeneous system in the thermodynamic Hmit [9^ . 
according to the well-known Hohenberg-Mermin- Wagner 
theorem. To avoid this technical difficulty, we therefore 
understand that the polarized gas under study is confined 
either in a box with a finite length L or in a harmonic 
trap (following the experiments) , although sometimes we 
would like to extend the length L to infinity. 

The key results of the present work are that the struc- 
ture of the ID FFLO state is clarified. The transition 
from the BCS state to the FFLO state is shown to be 
smooth, in marked contrast to the prediction of a first 
order transition in 3D [43|. Therefore, a ID phase sep- 
aration is excluded in the phase diagram of the uniform 
system. The phase separation in traps found in our pre- 
vious Letter is indeed simply an artifact of the paraboHc 
trap, as we anticipated. It is possible that similar effects 
are responsible for the phase separation observations in 



3 



the Rice experiment [2a|, which uses a high aspect ratio, 
elongated 3D trap. 

It should be emphasized that as well as being an in- 
structive theoretical test bed for the ground state prob- 
lem for a 3D g ID polarized Fermi gas in a trap can 
be realized exactly using two-dimensional optical lattices 
[9^ . [97I I . In these experiments the radial motion of atoms 
is frozen to zero-point oscillations due to a tight trans- 
verse confinement, while the axial motion is weakly con- 
fined. Thus, one can realize a low-dimensional quantum 
many-body system, and experimentally check the many- 
body predictions directly. This has also been recently 
carried out for a ID Bose gas [o^. lo^. 

The paper is organized as follows. In the following 
section, we outline the theoretical model for a ID spin- 
polarized Fermi gas. In Sec. Ill, we characterize the 
uniform phase diagram by using a simplified mean-field 
approach with a single-plane-wave like order parameter, 
i.e., the so-called FF solution for the FFLO state. This 
provides us with an approximate picture of the ground 
state of a ID polarized gas. An improved self-consistent 
BdG calculation is then given in Sec. IV, without any as- 
sumption for the order parameter. The underlying struc- 
ture of the FFLO states at all spin polarizations is then 
analyzed. The comparison between these two different 
mean-field approaches shows that the simple FF ansatz 
fails to capture the correct physics around the BCS- 
FFLO transition point. It therefore predicts the wrong 
type of transition. We conclude that in a 3D polarized 
gas case, the FF ansatz could lead to the same incorrect 
conclusion. In Sec. V the validity of these ID mean-field 
analyses in the weak-coupling or intermediate-coupling 
regime is checked using exact Bethe ansatz solutions. A 
quantitative phase diagram of a homogeneous gas is ob- 
tained by gathering all the information from these three 
methods. 

In Sees. VI and VII we study the trapped case, us- 
ing either the self-consistent BdG equations or the ex- 
act solution within the local density approximation. We 
again find a good agreement between these two results 
for weak and moderate couplings. The phase diagram 
of the trapped gas is thereby determined. We also cal- 
culate the local fermionic density of states of the FFLO 
states. A two-energy-gap structure is predicted, which is 
potentially useful for the experimental detection of FFLO 
states. Finally, Sec. VIII is devoted to the conclusions 
and some final remarks. 



II. MODELS 

Consider a polarized Fermi gas with a broad Feshbach 
resonance in a highly elongated trap formed using a two 
dimensional optical lattice [oHl- By suitably tuning the 
lattice depth, the anisotropy aspect ratio A = ujz/ujp of 
two harmonic frequencies can be extremely small. As 
long as the Fermi energy associated with the longitudinal 
motion of the atoms is much smaller than the energy level 



separation along the transverse direction, i.e., ^ 
hujp and NhuJz fii^p, where N is the total number of 
atoms, the transverse motion will be essentially frozen 
out. One ends up with a quasi-one dimensional system. 
The effective Hamiltonian of the ID polarized attractive 
Fermi gas the n niay b e described by a single channel 
model [2l|, M, [1Q3, [loH , 



+ Vtrap {x) - fJ-a (x) 



+ giD J dx^:l (x) (x) (x) (x) 



(2.1) 



where the pseudospins <J denote the two hyperfine 

states, and ^'o- (x) is the Fermi field operator that annihi- 
lates an atom at position x in the spin a state. The num- 
ber of atoms in each spin component is Na and the total 
number of atoms is N = + N^^. Two different chemi- 
cal potentials, = M ± <5/x, are introduced to take into 
account the population imbalance SN — N-^ ~ Ni > 0. 
The potential Vtrap {x) = defines a harmonic 

trap with an oscillation frequency lo — ujz va the axial 
direction. In such a quasi-one dime nsional geometry, it 
is shown by Bergeman et al. |l02l | that the scattering 
properties of the atoms can be well described using a 
contact potential giD^ix), where the ID effective cou- 
pling constant gm < may be expressed through the 
3D scattering length a^o, 



giD 



ma? [l- Aa^io/ap)' 



(2.2) 



Here Qp = ^JhJ(jrvjjp) is the characteristic oscillator 
length in the transverse axis, and the constant A — 
— C(l/2)/-\/2 ~ 1.0326 is responsible for the confine- 
ment induced Feshbach resonance |l02l . llQSl llQ4l | , which 
changes the scattering properties dramatically when the 
3D scattering length is comparable to the transverse 
oscillator length. It is also convenient to express g\B 
in terms of an effective ID scattering length, g^u — 
—2f? I (mai£)), where 



a\D 



p 



1-A 



> 0. 



(2.3) 



Note that in the definition of the ID scattering length, 
the sign convention is opposite to the 3D case. 

In this paper, we will assume a negative 3D scatter- 
ing length. In other words, the ID attractive polarized 
Fermi gas would be obtained experimentally from a 3D 
polarized gas onj^he BCS side of the Feshbach resonance 
magnetic field |il05i . 106]. 

In the absence of the harmonic trap, we measure the 
interactions by a dimensionless parameter 7, which is the 
ratio of the interacti on en ergy density eint to the kinetic 
energy density Ckm [l07l |. In the weak coupling limit, 
eint gion and Ckin ~ h?k'^/{2m) ~ h^ri^/m, where n 
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is the total linear density. Therefore, one finds 

mgiD _ 2 



7 = 



naiD 



(2.4) 



Thus, 7^1 corresponds to the weakly interacting limit, 
while the strong coupling regime is realized when 7^1. 

In the case of a trap, we may characterize the inter- 
actions using the dimensionless parameter at the trap 
center 70 = j{x — 0). For an ideal two-component Fermi 
gas with equal spin populations, the total linear density 



is 



Uideal (x) = no 1 - 



1/2 



in the Thomas-Fermi (TF) approximation, where 

2Ari/2 



no 

XTF 



(2.5) 

(2.6) 
(2.7) 



are respectively the center lin ear density and the TF ra- 
dius. Here aho = {muJz) is the characteristic oscilla- 
tor length in the axial direction. We thus estimate 



7o 



Ari/2 



a-ho 
aiD 



(2., 



In our previous Letter [85||, we have defined a dimen- 
sionless quantity k = Na\jj/a\^ to describe the interac- 
tions. These are related via 70 — T:/{^/n). 

Finally, we use a capital P = (N^ — Ni)/N to label the 
total spin polarization, and p — (n^ — ni)/n to denote 
the local (or uniform) spin polarization. 

To make the experimental relevance, we estimate the 
dimensionless interaction parameters for the on-going ex- 
periments on one-dimensional polarized Fermi gases. A 
gas of ^Li atoms in a three-dimensional optical lattice 
has been successfully produced by the MIT group [l9| . 
Thus, we consider the case of ^Li gas loaded into a two- 
dimensional optical lattice with the same parameters. 
Typically, in each one-dimensional tube the number of 
^Li atom is about N ~ 100. The transverse oscillator 
length Qp is related to the periodicity of the lattice d via 
Up = d/{ns^^'^) |l08l |. where s is the ratio of the lattice 
depth to the recoil energy. Taking s — A, the experi- 
mental value oi d = 532nm then yields Op ~ 120nm. An 
axial confinement of w ~ 2 7r x 400H z gives rise to an axial 
oscillator length aho = \J^I (moj) ~ 2/iTO. Further, the 
three-dimension scatt ering length of ^Li gas at the broad 
resonance is given by |l09l |. asd = — 1405ao[l -I- 300/ (B — 
834)] [1 + 0.0004(5 - 834)], where the magnetic field B 
is measured in Gauss and ao = 0.0529nm is the Bohr 
radius. We then use the relation. 



70 



1 



Ari/2 



at 



(1 - Aasn/ap) ' 



(2.9) 



to estimate the dimensionless coupling constant at the 
trap center. 
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Figure 1: (Color online) Dimensionless coupling constant at 
the trap center as a function of the magnetic field. This plot is 
designed specifically to represent a polarized gas of ®Li atoms 
in a two-dimensional optical lattice, assuming the same con- 
ditions as in the MIT experiment In detail, we take 
the total number of atoms as ~ 10^, and therefore in each 
tube the number of fermions is about A'^ ~ 100. The pe- 
riodicity of the lattice is d = 532nm, yielding a transverse 
scale Hp — 120nm. The axial confinement frequency ui ~ 
27rx400Hz, giving rise to an axial oscillator length a^o — 2fim. 
The 3D scattering length is related to the magnetic field via 
asd = -1405ao[H-300/(B-834)][l-h0.0004(B-834)], where 
the magnetic field B is measured in Gauss and ao = 0.0529nm 
is the Bohr radius. The dashed line in the figure shows the 
Feshbach resonance field. 



Fig. 1 gives the resulting 70 as a function of the mag- 
netic field B. We find that 70 ~ 0(1) above the Feshbach 
resonance. Throughout this work we shall take a cou- 
pling constant of 7 = 1.6. We note that there is already 
some indirect evidence for superfiuidity of a Fermi gas in 
a three-dimensional optical lattice at the magnetic 
field considered. On switching to a two-dimensional opti- 
cal lattice, the temperature in the experiments may still 
be low enough to generate the various one-dimensional 
superfiuid phases at zero temperature. 



Throughout the paper we shall mainly study two dif- 
ferent cases, either with a fixed total number of parti- 
cles and a fixed chemical potential difference, or with 
given numbers of both spin-up and spin-down particles. 
The system with two fixed chemical potentials may be 
considered as well. These three situations require the 
use of different canonical ensembles in thermodynamics. 
In the first two cases, we minimize the free energies of 
the system, Fsp{T, V, n, 6fi) and Fsn{T, V, n, 5n), respec- 
tively. While in the latter case, we minimize instead the 
thermodynamic potential fl{T, V, /j,, S^j,). 
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III. SINGLE PLANE WAVE APPROXIMATION 
IN A HOMOGENEOUS GAS 

We first consider a mean-field description with a single 
plane- wave FF type order parameter, to give the simplest 
qualitative picture of a homogeneous polarized Fermi gas 
[48]. At this point, we write the Hamiltonian l|2.ip in 
momentum space using a Fourier decomposition of the 
Fermi field operators. This results in: 



horn 



ka 



+91D C+ 2^j^^C+ 2_j^^Cp/2-fc'xCp/2+fc43.1) 

pkk' 

where — h^hF' /2ra is the kinetic energy. The single- 
plane- wave mean-field approximation amounts to decou- 
pling the interaction term using an order parameter 
A = -9iDY.k{(^q/2-kiCq/2+k-\) for the Cooper pairs, 
where we assume that the pairing occurs between a spin 
up atom with a momentum q/2 + k and a spin down 
atom with a momentum q/2 ~ k. As a result, the pairs 
possess a specific nonzero center-of-mass momentum q, 
whose value, together with the value of A, are to be de- 
termined. It is easy to see that after a Fourier trans- 
formation, the order parameter in real space acquires a 
single-plane- wave form, «.e., A(a:) — Aexp[iga;]. There- 
fore, within this approximation, we have a mean-field 
Hamiltonian, 



n 



MF 
horn 



giD 



Y (Cg/2-fciCg/2+fcT + 'l-C-) ■ (3-2) 

k 

Here, as a consequence of the constant linear density, 
Hartree terms like giDn~crcl„Cka merely introduce an 
overall shift for the chemical potentials. We indicate this 
by introducing the notation fla = ^J-a — gion-a for the 
shifted chemical potentials. 

To solve the mean-field Hamiltonian, it is conve- 
nient to use a Nambu spinor creation operator 'ip^ — 
(c^/2-i-feT' The Hamiltonian may then be rewrit- 

ten in a compact bilinear form, 

Whom = E M - A^) - Aa, + (e- - 5fi)] 4>k 

k 



5ii5"T"i + E (e/c - A + <5/i) , (3.3) 

giD k 



where 



(£9/2-1-*: ± eg/2-fe)/2, and and are the 
Pauli matrices. For convenience, we have defined. 



5jl — 5ii + 



gipn 
2 ' 
gmSn 



(3.4) 
(3.5) 



The bihnear Hamiltonian can be easily diagonahzed by 
working out the eigenvalues and eigenstates $^ of 



the two by two matrix [(e^ — fi)(Tz — Aax + (e^. 
Explicitly, we find that 



and 



^1 



Uk 
Vk 



6jl ± Ek, 



$1 



where Ek = [(4 - f^f + A^l^^^ and 



vl 



UkVk = 



1 - 

A 

2Ek 



4 - 

Ek 



Ek 



- 

(3.6) 
(3.7) 

(3.8) 
(3.9) 
(3.10) 



From the eigenstates <i>J, it is natural to define Bogoli- 
ubov quasiparticle operators, which are given by: 



Uk, 



ipk 



The bilinear mean-field Hamiltonian then becomes 

k 

'Y.\^k 



(3.11) 



njMF 

'^hom 



gion^ni + ^{e+ - jl- E^ 
7 ~5fi\ a+^ak] 

k 



(3.12) 



The thermodynamic potential is obtained by replacing 
a^^afeo- by its thermal statistical average values, i.e., the 
Fermi distribution function f{E^) = l/(exp[/3i?^] 4- 1) 
with (3 = l/{kBT) as the inverse temperature. At zero 
temperature where (3 goes to infinity, the Fermi distri- 
bution function f{x) reduces to a step function Q (—a;), 
i.e., 9 (x > 0) — 1 and 9 (x < 0) = 0, so the resulting 
thermodynamic potential has the form: 



A2__ 

giD 



gion^ni + ^ 

A; 

t + — <^m] 9 I 



/i 



Eh 



-El 



+ T.[^^-'k+5A^{-Ek)^ (3-13) 



The values of the order parameter A and of the pair- 
ing momentum q are determined by finding the station- 
ary points in the (A, q) plane of the thermodynamic po- 
tential, i.e., dil/dA — and dVl/dq — 0, with given 
chemical potential difference (S/i, or the requirement of 
number conservation, 5n = —d^l/dS^j,. This gives us two 
distinct procedures for defining the mean-field solution. 
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Figure 2: (Color online) Landscape of the zero-temperature 
thermodynamic potential of a uniform gas at an interaction 
strength 7 = 1.6. Here, we take a single-plane-wave approx- 
imation for the order parameter, and normalize it using the 
full gap of an unpolarized Fermi gas, Ao = 0.34658eF, where 
ep is the Fermi energy. The chemical potential is fixed at 
jl = 1.04594eF. The competing ground states are (i) a nor- 
mal Fermi gas with A = 0, (ii) a fully paired BCS superfluid 
with A = Ao, q = 0, and 5n = 0, (iii) a finite momentum 
paired FF superfluid with A < Ao, g / 0, and Sn 7^ 0, (iv) 
a breached pairing or Sarma superfluid with A < Ao, q — 0, 
and Sn 7^ 0, and (v) a saddle point phase intervening between 
the local BCS and FF minima. We note that the last two 
phases are unstable with respect to phase separation. 



analogous to the grand-canonical (fixed chemical poten- 
tial difference) and canonical (fixed number difference) 
ensembles in thermodynamics. 

Once these variational variables are obtained, we cal- 
culate straightforwardly the total free energies = 
il+jin = Fsi_i.+giD{'n^ +Sn^)/'i or Fgn — fl+^n+S^Sn — 
Fsn + gioin^ — 5ii?)/A of the gas, depending on whether 
the chemical potential difference (5/i or the number differ- 
ence 5n = ni — ni is fixed, as indicated by the subscript. 
Note that at zero temperature the value of the free energy 
Fsn is equal to the total ground state energy E. We have 
also defined two free energies Fg^ and Fgn in the absence 
of the Hartree terms. In the detailed calculations, for 
a uniform system we take respectively the Fermi energy 
ep = h?kp/{2m) and the Fermi wave vector kp — Trn/2 
(of a unpolarized ideal gas ) as the units of the energy 
and of the momentum, by letting h = 1 and 2m = 1. 



A. Qualitative phase diagrams 

Generally, there are several possible stationary solu- 
tions in the landscape of the thermodynamic potential. 
On the weak coupling side we find only three stable com- 
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Figure 3: (Color online) Comparison of the free energies 
of Fsfi available mean-fleld solutions at a coupling constant 
7 = 1.6 and at zero temperature, with the free energy of the 
normal gas Fjv being subtracted. With increasing the chemi- 
cal potential difference, the gas turns from a BCS superfluid 
to a FF superfluid at 5jl ~ O.68A0, and flnally becomes a 
normal gas above Sp. = 2eF. 



peting ground states, corresponding to local minima of 
the landscape. As shown in Fig. 2 for a coupHng con- 
stant 7 = 1.6, these are the unpolarized (BCS), par- 
tially polarized (FF), and a fully polarized or normal 
(N) phases. The other two states, denoted as "Sarma" 
and "saddle point" phases, are unstable with respect to 
phase separation [43|. Note that in the figure, the order 
parameter A and the center-of-mass momentum q are 
measured in units of the full gap of an unpolarized gas, 
Ao ~ 0.34658e_F. We have fixed the chemical potential 
at its unpolarized value, p. ~ 1.04594ei?, and have taken 
the chemical potential difference to be 5jl = O.TSAq. 

For an interaction strength 7 = 1.6, the evolution of 
the ground states with increasing chemical potential dif- 
ference is given in Fig. 3. Here we search for the ground 
state by minimizing the free energy . As Sjl increases 
from zero, the free energy of the BCS state is initially 
lowest, but rises very rapidly. It intersects with that of 
the FF state at about dp, = O.68A0. A first order quan- 
tum phase transition then occurs in mean-field theory, 
since the first order derivative of free energies at the in- 
tersection point is discontinuous. The apparent hystere- 
sis (presence of the FF state before the transition point) 
is also the mark of a first order phase transition. After 
that, the free energy increases slowly towards the nor- 
mal state value. Precisely at Sp = 2e_F, the gas enters 
smoothly into a fully polarized normal state, where the 
spin polarization p — (n| — ni)/{n^ is strictly equal 

to one. Hence, differing from the 3D situation, a partially 
polarized normal phase is excluded in ID. We present, re- 
spectively, the value of the order parameter and the spin 
polarization as a function of the chemical potential dif- 
ference in Figs. 4a and 4b. The first order transition 
from BCS to FF states becomes much apparent due to 
the jump of the order parameter and of the spin polar- 
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Figure 4: (Color online) Evolution of the mean-field (FF) 
order parameter and of the spin polarization, with increasing 
chemical potential difference. The arrows point to the phase 
transition positions. The parameters are the same as in Fig. 
3. 



ization. We will show later, however, that this apparent 
first order transition is simply an artifact of the single- 
plane- wave approximation for the order parameter. 

By changing the coupling constant, we can determine 
a phase diagram in the plane of the interaction strength 
7 and chemical potential difference Sfi, as shown in Fig. 
5a. The solid and dashed lines separate the FF state 
from the normal and BCS phases respectively, and con- 
verge to a single curve above 7 ~ 7. Converting the 
chemical potential difference to a number difference, we 
obtain a phase diagram in the 7 — p plane in Fig. 5b. 
The area under the dashed line has no correspondence in 
Fig. 5a and belongs to the "saddle point" solution, which 
is unstable towards phase separation. This may be the 
precursor of a phase separation phase. Overall, all the 
basic features found here are qualitatively similar to that 
in 3D El. 



B. Analytic results in limiting cases 

We discuss some analytic results that can be obtained 
in the weakly interacting limit of 7 0. The simplest 
one is the unpolarized BCS state, for which the chemi- 
cal potential /2 is essentially the Fermi energy ep. The 
stationary condition dfl/dA = then leads to a gap 



1.0 



0.8 



0.6 



0.4 



' 1 

(a) 


' 1 ' 1 






FF 


\ ^ 
\ 


^ ~ > _ 


BCS 

1 1 


1,1,1, 


1 


2 3 


4 5 6 7 


' 1 

(b) 


' 1 ' 1 

N 


1 ' 1 ' 1 ' 




FF 




1 


2 3 


4 5 6 7 



1.5 



1.0 



0.5 



0.0 



Figure 5: (Color online) (a) Phase diagram in the plane of 
the interaction strength and the chemical potential difference. 
Within the single-plane- wave assumption for the order param- 
eter, the transition from a BCS superfiuid to a FF state is of 
first order (dashed line), while from a FF state to the normal 
state it is continuous (solid line), (b) Interaction strength vs 
polarization phase diagram. The shadow region is unknown, 
and presumably is an artifact of the single-plane- wave approx- 
imation. 



equation, 
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(3.14) 



The integration can be worked out analytically for small 
Aq. One finds that 



Sep exp 



(3.15) 



analogous to the standard 3D BCS result Aq^ ~ 
Sep e'xp[Tr/ {2k pa) — 2]. For the FF state at a large chem- 
ical potential difference, the value of the order parameter 
is even smaller. To a good approximation, we find that 



qkp 



epA 
8jl, 



Aep '■ 



and hence: 
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(3.16) 
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(3.18) 
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intriguing two-energy-gap structure. The midgap state 
around e = is a salient feature of the spatially modu- 
lated order parameter jo^]- 



IV. SELF-CONSISTENT BDG IN A 
HOMOGENEOUS GAS 

We now turn to a more realistic mean-field calcula- 
tion without resorting any approximation for the form of 
the order parameter. We consid er th e BdG equations of 
the ID polarized Fermi gas [Hll, \lld^ . starting from the 
Heisenberg equation of motion of the Hamiltonian (|2.ip 
for {x,t) and i {x,t) (without the trap potential): 



Figure 6: (Color online) Local fermionic density of states of 
a uniform polarized Fermi gas, with a single-plane-wave form 
for the order parameter. Note that there is a prominent two- 
energy-gap structure in the FF state. 



From the prefactor, the order parameter A vanishes ex- 
actly aX Sp, = 2eF- At the same time fi — 2eF, indicating 
that the FF state changes smoothly into a fully polarized 
normal state. 



C. Local fermionic density of states 

The Bogoliubov quasiparticle amplitudes {uk, Vk) and 
energy Ek appear in the zero temperature spectrum of 
the single fermionic excitations. We characterize the ex- 
citation spectrum using the local fermionic density of 
states, Per (e), given by 

PT(e) = J2^lH^-K) + Y.''lH^-Ek){s-i^) 

k k 

p,{e) = J2vlS{e + E+)+Y,ul6{e + E^) (3.20) 

k k 

For an ideal gas with equal populations, the density of 
states can be calculated analytically. 



bk / \ bk / \ 

P] (e) = Pi (e) 



V2rn 1 
2'Kh y/e + p. 



(3.21) 



which we have regarded as a background density of states. 
It has a band edge (square root) singularity at e = —p,. 

We plot in Fig. 6 the local density of states for a 
one-dimensional BCS superfluid, and the FF phase at 
p = 0.12, as well as the background density of states. In 
an FF state, the spin up and down density of states are 
exactly the same, but are shifted downwards or upwards 
respectively by an amount Sp. For clarity, in the figure we 
show only one branch, i.e., the spin up density of states 
after an upwards shift. Compared to the BCS superfluid, 
the local density of states of the FF phase exhibits an 



ih- 
ih- 







dt 


2m 






dt 


2m 



Mi 



*T + .9ir'*t*i*T' (4-1) 
*i - .9iD4'+4'i*T- (4-2) 



Within the mean-field approximation, we replace the 
terms gio'^^'^ i^] and gm^l^^ i^] by their respective 
mean-field decoupling 



^ -A(a;)*+ + giDni(:x)^^, (4.3) 



and 



-A(a;)^'+ +5-iDnt(x)^'x, 



(4.4) 



where we have defined an order parame- 
ter A(a;) — — gi£)(5'j^(x)4'|(a;)) and densities 
n^{x) = [x)^ cr{x)) . The above decoupling 
thus yields. 



ih 



ih 



(4.5) 



dt 



[?^J-Mi]*i + A(x)*+, (4.6) 



where H% = —h^V'^/ {2m) + giong- (x). We solve the 
equation of motion by inserting the standard Bogoliubov 
transformation: 

*T = EKt (x) c.Te-^^"^*/'' + v;^ (x) c+ e^^"^*/'-^, 
n 

= EKi(^)<i^^''''^*^'-«'^T(^)c.Te-^''''^*%-7) 

This gives rise to the wel l-kno wn BdG equations for the 
Bogoliubov quasiparticle 



- p, -A{x) 
-A*{x) -n% + ps 



^7] (J 



7](T 



where the wave functions m^o- {x) and Vria {x) are normal- 
ized by 



dx \Una {x)f' + \Vj^a {x)\^ =1, (4.9) 
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and EritT is the corresponding excitation energy. 

We note that the unequal chemical potentials of spin 
states in the BdG equations break the particle-hole sym- 
metry. This leads to different quasiparticle wave func- 
tions for the two components. However, one may easily 
identify a one to one correspondence between the solution 
for the spin up and spin down energy levels, i.e., 



E„ 



-En 



and 







-~V*^a {x) 


Vrja {x) 







(4.10) 



(4.11) 



Because of this symmetry of the BdG equations, there- 
fore, we may consider the spin up part only. Letting 
Ur, (x) = (x) and w,, (r) = (x), we then remove 
the spin index in the equations. 



are two branch solutions for the quasiparticle energy 

= {<^q/2+k-<^q/2-k)/'2.-5jl+Ek Snd E^ = {tq/2+k- 

£g/2-fc)/2 — ^/i— with the corresponding quasiparticle 
wave functions, 



Vk 



Uk 
Vk 



and 



Uk 
Vk 



Ek — E, 



-Vk 



^k ' 



(4.21) 



(4.22) 



respectively, exactly the same as in Eqs. I|3.6p and p.7p . 
Accordingly, the linear densities take the form, 

(x) = <f ^^t) + E ^^k) . (4-23) 

k k 

niix) = J2vlf {-E+) +J2ulf {-E-) ,(4.24) 



m 



-A*ix) -ni + fii 



Un [x) 

Vr, (x) 



En 



Urj (x) 
Vrj {x) 



(4.12) 

The order parameter A(a;) and the linear number den- 
sities Ua- (x) should be determined self-consistently, ac- 
cording to their definitions, respectively. 



n^ix) = J2<(''>vi^)fiEv), (4-13) 
'^i(^) = E"'5(^)^')(^)-^("^'))' (4.14) 

A{x) = -g,DY.M^Kix)fiEv)- (4.15) 
»? 

where the summation runs over all the energy levels, in- 
cluding these with negative energies E^ <0. 

We note also that the single-plane-wave approxima- 
tion described in the last section can be recovered by 
replacing the level index "77" with a wave vector k, and 
approximating. 



Unix) 

Vjjix) 



En 



Uk exp 
Vk exp 
Ek, 



-^[^ + k) X 
-z(|-fc) x_ 



(4.16) 

(4.17) 
(4.18) 



so that the order parameter reduces to 



A(a;) = -giDj2ukVkf{Ek)eyi-p[iqx] = Ac-Kp[iqx], 

(4.19) 

and the BdG equations become. 



-A 



-A 

-^q/2~k + Mi 





Uk 


= Ek 


Uk 


For the 











(4.20) 

where as before, we have used the notations /i| = 
MT ~ dioni and Mi = A*i ~ dionp Apparently, there 



which turn out to be position independent due to the 
plane-wave form of the wave functions. 



A. Hybrid BdG strategy 

We apply the above BdG formalism to a uniform Fermi 
gas with finite atoms. To this end, we consider a gas of 
N fermions in a box of length L using periodic boundary 
conditions, i.e., the underlying wavefunction (p{x) satis- 
fies ip{x — +L/2) — ip{x — —L/2). The small boundary 
effect due to the finite size of L could be weakened or 
removed by enlarging the value of L. 

In any practical calculation, because of the computa- 
tional limitations, the summation over the quasiparticle 
energy levels in Eqs. (|4.13p . I|4.14p and l|4.15p must be 
trun cate d. We therefore following the idea of Reidl et 
al. |lll| develop a hybrid approach with the introduc- 
tion of a high-energy cut-off E^, below which we solve 
the discrete BdG equations. Above the cut-off, we use a 
semiclassical plane- wave approximation for the wavefunc- 
tions, which should work well for sufficiently high- lying 
states. 

The first step toward solving the discrete BdG equa- 
tions is to assume a real order parameter A (a;) and then 
expand the quasiparticle wavefunctions u {x) and v (x) 
using a complete basis of single particle wavefunctions in 
the box ifnix) with energy levels e„ (n = 0, 1, 2, ...), i.e., 



u{x) = A„y?„(a 

n 

v{x) = ^B„(p„(a 



(4.25) 
(4.26) 



fn{x) 



■\/2/i cos [mrx/L] , if n is even; 
/2/Lsin [{n + 1) irx/L] , if n is odd; 



(4.27) 
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and 

ft^TT^n^/ (2rnL^) , if n is even; 
h^n^ {n + if / (2mL2) , if n is odd; ' ^ 

The solution of the BdG equations then becomes a matrix 
diagonalization problem, 



where Ek{x) = y/ (ck — A*)^ + i^) and the subscript 
"c" means the continuous contribution from high-energy 
levels. 

The discrete and continuous parts of the order param- 
eter may be combined together to give, 



'TVOT _j_ A/Y I 

nn' ' nn' 



-Ann' 



Mi 





An' 


^ E 


An 




Bn' _ 


Bn 



V/Ji^) E M^Ki^)fiEv), (4.42) 

|B„|<-Ee 



where the matrix elements 
K 



^nn' 



(4.29) 
(4.30) 



where we have defined a position dependent effective ID 
coupling constant g^^J [x), which satisfies. 



(^n Mcr) ^nn' : 

+L/2 

Mnn' = 9lD J dxifn{x)nif {x) ipn>{x), (4.31) 
-L/2 



+L/2 

Ann' = / dxipn{x)A{x) (p„'{x). 
-L/2 



(4.32) 



The coefficients of the eigenstate has to satisfy the 
condition J^n i^n + B^) — 1 due to the nor- 
malization of the quasiparticle wavefunctions, i.e., 
Cl'^^dx {u\x)^v\x)\^\. 

These discrete spectra (labeled by an index "ry") con- 
tribute to the linear densities and the order parameter as 
follows, 

n^i{x) = u;{x)ur,ix)fiE^), (4.33) 

I I < Ec 

n^,{^) = E v;{x)vr,ix)f{-E,), (4.34) 
Ad{x) = -giD E '^v(^)'"vi^)fiEv)^ (4-35) 

where the subscript "d" refers to the discrete levels. 

On the other hand, for the h igh-l ying states we take 
the semiclassical approximation |lllj . 

Uri{x) u{k, x) exp[ikx] , (4.36) 
Vjjix) v{k,x) e:xjp[ikx] , (4.37) 
^ E{k), (4.38) 

where we have regarded the wavefunctions locally at po- 
sition X as plane waves, whose amplitudes u{k,x) and 
v{k, x) are normalized according to u^(fc, x) + w^(fc, x) — 
1. Keeping the most important pair correlation terms 
only, it is straightforward to show that at low tempera- 
tures, 



n|c [x) 
nic (x) 
Ac {x) 



k 



2Ek{x) 

<^k- ^J■ 

2Ek{x) 
A(x) 
't 2Ek{x) 



e[Ek{x) + Sfi- Ei^,39) 
e[Ekix)-S^i-EQLAO) 
e[Ek{x)+Sti-Ec](,4.41) 



1 



1 



91d (x) giD 



where 



(4.43) 



(4.44) 



The summation over the momentum k may be converted 
into a continuous integral of the energy. As a result, we 
obtain. 



n-fc (x) 



(2m) 



1/2 



de 



^{e-Siif-A-^ {x) 



nic (x) 



fi+^{e-Sf,f-A^ (x) 
(2m)^/^ 



1/2 ■ 



(4.45) 



de 



^{e + Sf,f-A^ (x) 



fi+^Jie + Sfif^A^ (x) 



1/2 ■ 



(4.46) 



and 



(2m) 



1/2 



de 



^{e-S^,f^A^ (x) 



^i+Jie-Sf^r-A^ (x) 



1/2- 



(4.47) 



We can now summarize the entire procedure used to 
obtain the BdG solutions. The key step is to solve the 
eigenvalue problem (|4.29p . As the calculation of matrix 
elements involves the order parameter and linear densi- 
ties that are yet to be determined, a self-consistent itera- 
tive procedure is required. For a given number of atoms 
{N ^ N-^ + Ni and SN ^ - Ni), temperature and 
interaction coupling gm, we: 

(a) start with an initial guess or a previously determined 
better estimate for A (x), 



11 



(b) solve Eqs. (|4.43p and l|4.47p for the effective coupling 
constant, 



(c) then solve Eq. I|4.29p for all the quasiparticle wave- 
functions up to the chosen energy cut-off to find 
Uri{x) and Vri{x), and finally determine an im- 
proved value for the order parameter from Eq. 

(031. 

During the iteration, the density profiles n-^ix) — 
n]-dix) + n'^c{x) and ni{x) — nid{x)+nic{x) are updated. 
The chemical potentials /i and Sfi are also adjusted 
slightly in each iterative step to enforce the number- 



conservation condition that /I^^/^^ dx[n^{x) + ni{x)] 



--N 



and J^^l^ dx[n-^ (x) 
gence is reached. 



{x)]—SN, until final conver- 



B. The structure of FFLO states 

Using the self-consistent BdG formalism we can work 
out the detailed structure of mean-field or FFLO states. 
To make the equations dimensionless, as before we take 
the Fermi wave vector kp = Trn/2 — TrN/{2L) and the 
Fermi energy ep = h'^kp/{2m) as the units of the mo- 
mentum and energy, respectively, i.e., by setting h = 1 
and 2m = 1, and kp = 1. Therefore, the size of the 
box L — 'kN/2 can be enlarged by increasing the number 
of total atoms TV. In the following calculations, we use 
N — 200, which in most cases we find is large enough 
to effectively minimize the boundary effects. Further, we 
take a cut-off energy Ec — 16ep. This cut-off energy is 
already sufficient large because of the high efficiency of 
our hybrid strategy. Accordingly, we setup a set of single- 
particle-state basis ipn{x), with the highest energy level 
larger than the cut-off energy. 

The initial guess for the order parameter A(a;) could 
be arbitrary. However, we find that in general there 
are many locally metastable solutions after the iteration, 
which can be classified uniquely by their periodicity. This 
is due to the existence of the periodic boundary condition 
that requires that the order parameter should be a peri- 
odic function of length L/n, where n is an integer. We 
therefore compare the energy (or free energy) of the so- 
lutions with different periodicity, and select the one with 
the lowest energy as the ground state. 

We present in Fig. 7 the spatial distribution of the 
order parameter A (x) and the local spin polarization 



P{x) 



(4.48) 



for a uniform Fermi gas with total polarization p = 0.03 
(a) and p = 0.16 (b) at a typical coupling constant 
7 — 1.6. The most notable feature of the figure is that 
at a small total polarization (Fig. 7a), the order pa- 
rameter switches between two values: -I-Aq and — Aq, 
where Aq is the full gap of an unpolarized gas at the 
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Figure 7: (Color online) Spatial structures of the inhomoge- 
neous FFLO states at an interaction strength 7 = 1.6 and 
at two spin polarizations as indicated. The calculations have 
been done for a uniform gas confined in a box, using the self- 
consistent BdG equations. The solid line and the dashed line 
refer to the order parameter and the local spin polarization, 
respectively. 



same coupling. Many instantons and anti-instantons (or 
kinks and anti-kinks) then appear and carry the excess 
spin up (majority) atoms since the local polarization p{x) 
shows pronounced peaks right at the position where the 
order parameter vanishes. These features are not unlike 
a phase separation, except that a regular, periodic do- 
main structure is obtained. Thus, in the limit of small 
polarization, the order parameter may be viewed as an 
instanton gas, with the number of instantons roughly pro- 
portional to the spin polarization. Within this picture, 
we anticipate that an FFLO state emerges as soon as the 
polarization becomes nonzero. In contrast, for a large 
total polarization (Fig. 7b), the order parameter is well 
approximated by a cosine function, as expected earlier 
by Larkin and Ovchinnikov. It is a superposition of two 
single-plane-waves going in opposite directions, with a 
much reduced amplitude compared to Aq. 

We note that in the weak coupling limit, a snoidal 
solution of the order parameter for the BdG equations 
was found analytically if one linearizes the single parti- 
cle spectrum at the Fermi surface [H, |9^, which gives 
qualitatively the same behavior as shown in Fig. 7. 
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Figure 8: (Color online) Spin polarization versus the chemical 
potential difference at an interaction strength 7 — 1.6, ob- 
tained from the single-plane-wave approximation (solid line) 
and the self-consistent BdG calculations (open circles). While 
the spin polarization in the FF state shows a jump as a func- 
tion of the chemical potential difference, the more accurate 
self-consistent BdG prediction suggests that the spin polar- 
ization emerges from zero continuously with increasing the 
chemical potential difference. The dashed line is a power-law 
fit to the self-consistent BdG results. 



C. Phase diagram from BdG solutions 

We examine the phase diagram obtained by the single- 
plane-wave approximation (Fig. 5). For this purpose, we 
compare the results of the spin polarization versus the 
chemical potential difference, as predicted respectively by 
the self-consistent BdG formalism and the single-plane- 
wave approximation or the FF solution. As shown in Fig. 
8, the self-consistent prediction agrees very well with that 
of the FF solution at a large chemical potential differ- 
ence. However, approaching to the BCS-FFLO transi- 
tion point, they differ largely. The quick fall of the spin 
polarization in the self-consistent BdG indicates strongly 
the existence of a FFLO state with an arbitrary small 
spin polarization. As the spin polarization is a first or- 
der derivative of the energy, this is a solid evidence for 
the smooth transition from the BCS state to the FFLO 
state. We therefore conclude that although the single- 
plane- wave approximation gives a reasonable description 
at the large chemical potential difference, it does not pre- 
dict the correct phase transition between BCS and FFLO 
states. 

We may extract the critical behavior at the transition 
point by numerically analyzing the self-consistent data. 
Assuming a pow-law dependence of the spin polarization 
on the chemical potential difference, p oc {Sfj, — Sfic)", 
we find that a ~ 0.4, in good agreement with a non- 
perturbative bosonization prediction j9^, a = 0.5. The 
small discrepancy may be caused by the use of a finite 
length L, which becomes increasingly in-efficient due to 
the divergent correlation length towards the transition 
point. 



Figure 9: (Color online) Local fermionic density of states of 
a uniform polarized Fermi gas at an interaction strength 7 = 
1.6, calculated using the self-consistent BdG equations. 



D. Local fermionic density of states 

We finally calculate the local density of states in the 
self-consistent BdG solutions, which is given by, 

= Y.ul{x)S{e~Er,), (4.49) 
n 

Pi(x,e) = J2''v(=^)H^ + E^)- (4.50) 

n 

In Fig. 9, we show how the local density of states at origin 
evolves with increasing the spin polarization p from zero 
to 0.12. Here a small spectral broadening of about 0.02e_F 
has been used to regularize the delta function. We find 
again a nonzero density of states at the Fermi surface for 
a polarized Fermi gas, contributed by the mid-gap states. 
As a result, the original BCS gap of a width 2Ao is split 
into two sub-gaps with a much smaller width. 



V. EXACT BETHE ANSATZ SOLUTION IN A 
HOMOGENEOUS GAS 

The validity of mean-field results in ID is not immedi- 
ately clear, as pair fiuctuations become increasingly im- 
portant in lower dimensions. Fortunately, without the 
trap the Hamiltonian (|2.ip of a free polarized Ferr ai g as is 
exactly soluble, using the Bethe ansatz technique (sgLIstI. 
We therefore can use the exact solution as a benchmark 
to test the validity of various mean-field approaches. 

In the thermodynamic limit, the ground state of a ho- 
mogeneous gas with fixed linear densities and may 
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be obtained from a set of Gaudin integral equations |87l |. 
1 B 
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where e^s is the ground state energy density, the cou- 
pUngs c = nj and c' = c/2. The functions p{k) and cr(A) 
are, respectively, the quasi-momentum distributions with 
the cut-off rapidities Q and B to be determined by the 
normalization condition for Sn — n-^ — ni and n|. The 
last term in egg is simply the contribution from rij paired 
two-fermion bound states with binding energy 
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ma 
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ID 



The chemical potential and the chemical potential differ- 
ence can be obtained by /i = dcgs/dn and Sp — degs/dSn, 
respectively. 



A. Gaudin solutions 

The Gaudin integral equations have to be solved nu- 
merically for a general spin polarization p — 5n/n. To 
do so, we introduce two new variables x ~ k/Q and 
y = A/B, and rewrite the quasi-momentum distribution 
functions. 
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9s (y) 
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Further, the two cut-off rapidities may be represented 
by, respectively, Q = wf/Xc and B = nj/Xs- In such a 
way, the Gaudin integral equations can be rewritten in a 
dimensionless form. 
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together with the normalization conditions, 
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Figure 10: (Color online) Gaudin solution for the dimen- 
sionless quasi-momentum distributions at a spin polarization 
p = 0.5 and at several interaction couplings as indicated. 



Numerically, the dimensionless integral equations have 
been solved by decomposing the integrals on a grid with 
N points {xi;xi G [— l,-)-l]}. In detail, we start from 
a set of trial distributions g''c\xi) and gf\xi)^ and the 
corresponding parameters of A^"^ a nd a I"^ . Following the 
standard method for the integrals |l07l |. we obtain gdxi) 

(0), 



and gsixi). Let gc^\xi) = ag^^^Xi) -I- (1 - a)gcixi) and 
gi-^\xi) — agf\xi) -I- (1 — a)gs{xi) (where a is a posi- 
tive real number between and 1, depending the value 
of the spin polarization) be the new trial distributions, 
and update x[^^ and Ai^"* accordingly. Repeat the above 
procedure until gdxi) and gs{xi) agree with their trial 
distributions within a certain range. Then, the energy 
density 



2^3 



"2^ 



■e(7,p) -n^eb 



(5.12) 



is calculated by: 



^3 +1 3 -1-1 

e (7,P) T3 / x'^gc {x) dx + ^ J 2x^gs (x) dx. (5.13) 

We find that this iterative method for solving the Gaudin 
integral equations is very stable. The chemical potential 
and chemical potential difference can also be calculated 
accurately by a numerical derivative. 

For an illustrative purpose, we plot in Fig. 10 the 
quasi-momentum distribution functions gs{x) (Fig. 10a) 
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Figure 11: (Color online) Gaudin solution for the dimension- 
less parameters Ac and As, as a function of the spin polariza- 
tion at an interaction strength 7 — 1.6. 



Recalling that rij = n(l — p)/2, the chemical potential, 
as well as the first two terms on the right-hand side of 
the chemical potential difference, coincide in magnitude 
with the chemical potenti al of a Tonks- Girardeau bosonic 
gas of paired nj^ dimers [l07l |. which is fermionized due 
to strong attractions. The third term in the chemical 
potential difference, on the other hand, is equal to the 
chemical potential of residual unpaired n-\ — fermions. 
Therefore, in the strong coupling regime the polarized 
gas behaves like an incoherent mixture of a molecular 
Bose gas and a fully polarized single-species Fermi gas. 

The analytic derivation in the weak coupling limit 
7 <C 1 is much more subtle since the quasi-momentum 
distribution gdx) contains a sharp jump whose width 
7) is extremely small, as shown in Fig. 10b for 
7 = 0.016. However, as a leading approximation, we 
may take gc{x) as a step function. It is then easy to 
show that (7 ^ max{p, 1 — p}). 



and gc{x) (Fig. 10b) at a spin polarization p — 0.5 for 
three interaction strengths as indicated. As 5s (x) and 
gdx) are both even functions, we show only the part with 
a positive x. For a large interaction strength, they ap- 
proach I/tt and l/(27r) respectively. On the other hand, 
for a weak interaction, gs{x) reduces to l/(27r) and gc(x) 
jumps from zero to l/(27r) at a certain value of x. 

At 7 = 1.6 the dimensionless parameters Ac and As as 
a function of the spin polarization are shown in Fig. 11. 
They diverge respectively as 1/p and 1/(1 — p) when the 
spin polarization goes to or 1. 



B. Analytic results in limiting cases 

The asymptotic behavior of the Gaudin solution may 
be obtained in the strongly and weakly interacting limits. 
For a strongly interacting gas, for which the dimension- 
less coupling constant 7^1, the parameters Ac and As 
are sufficient large. Therefore, the integrals in the Gaudin 
equations becomes extremely small. Hence, the quasi- 
momentum distributions gc(x) and gs{x) are essentially 
constant. Expanding to the order 1/7"^, we find that. 
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It is then straightforward to show that to leading order 
in 1/7, 
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As a result, the ground state energy density and the 
chemical potentials are given by 



e{l,p) 
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(5.21) 
(5.22) 
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where the first term on the right-hand side corresponds 
to an ideal polarized gas, while the second term arises 
from the mean-field Hartree-Fock interactions. We note 
that a non-perturbative term of order 7^ In 7 will occur if 
one improves the quasi-momentum distribution functions 
by explicitly taking into account the width of the jump 
in gc{x). 



C. Mean-field approaches versus exact solutions 

We are now ready to verify the accuracy of the mean- 
field approaches. In Figs. 12 and 13, we compare the 
energy and chemical potentials of the exact Gaudin so- 
lutions with that from mean-field calculations, with ei- 
ther a single-plane-wave like (labeled as "FF") or a self- 
consistently determined (denoted by "S'C-BdG") order 
parameter. For comparison, the energy of an ideal po- 
larization gas is also shown. For a moderate interaction 
coupling 7 = 1.6, we find a reasonable agreement. The 
residual discrepancy could be ascribed to pair fiuctua- 
tions, which are small but not neghgible. We have also 
checked that the agreement becomes increasingly bet- 
ter (as expected), with decreasing interaction strength. 
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Figure 12: (Color online) Comparison of the mean-field en- 
ergy to the exact results obtained from the Bethe ansatz so- 
lution at an interaction strength 7 — 1.6. For a reference, we 
plot also the energy of an ideal polarized gas. Presumably, the 
small discrepancy between the mean-field and exact results is 
due to the pair fiuctuation effects. 




Figure 14: (Color online) Phase diagram of a one-dimensional 
homogeneous spin-polarized Fermi gas. The dot-dashed line 
refers to the asymptotic expression of the critical chemi- 
cal potential difference in the weak coupling limit, i.e., Eq. 
I|5.27p . while the two dashed lines are respectively, the strong- 
coupling expansion of the critical chemical potential differ- 
ence, as described in Eqs. I|5.24p and l|5.25p . 
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Figure 13: (Color online) Comparison of the mean-field chem- 
ical potentials to the exact results obtained from the Bethe 
ansatz solution at an interaction strength 7 — 1.6. The arrows 
point to two critical chemical potential differences, between 
which a polarized superfiuid exists. 



With these observations, we therefore confirm the vahd- 
ity of the mean-field theories for the weakly and moder- 
ately interacting regimes. 

On the other hand, the good agreement between the 
Gaudin solutions and the mean-field results suggests 
strongly that the partially polarized solution found in 
the exact Bethe ansatz method is of FFLO character. 
We note that a calculation of the nonlocal pair correla- 
tion functions in the exact solution would be very useful 
to unambiguously determine its structure. However, this 
is extremely difficult due to the complicated ground state 
wavefunctions from the Bethe ansatz. 



D. Quantitative phase diagram of a homogeneous 
polarized Fermi gas 

Gathering all the information from the Gaudin integral 
solutions and the two mean-field results, we arrive at a 
quantitative phase diagram for a homogeneous polarized 
Fermi gas [ssl. [93|. For a given interaction strength the 
chemical potential difference takes values between two 
thresholds, S^c,p=o and 5^c,p=i, as indicated by arrows 
in Fig. 13 for 7 = 1.6. Below the first threshold Sfic,p=o, 
the gas persists in the BCS-like superfiuid state with zero 
polarization (SF), while above the second critical value 
<5/ic,p=i, a fully polarized normal state appears (N). In be- 
tween, a superfiuid state with finite polarization (SFp) 
is favored. As stated earlier, the SFp has a FFLO struc- 
ture in character. Physically 5^c.,p=o is the energy cost 
required to break spin-singlet pairs in unpolarized super- 
fiuid, i.e., the spin gap, while Sfic,p=i is also associated 
with the pair-breaking (for the last pair) , but is enhanced 
due to the Pauli repulsion from existing fermions. The 
dependence of Sfj,c,p=o and S^c,p=i on the parameter 7 is 
reported in Fig. 14, constituting a homogeneous phase 
diagram. 

The behavior of the critical chemical potential dif- 
ference in the weak and strong couphng limits may be 
worked out analytically. In the strongly interacting 
regime of 7 ^ 00, from its asymptotic expression (|5.18p 
we find that. 
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While in the weakly interacting limit of 7 ^ 0, only 
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Figure 15: (Color online) Same phase diagram as in Fig. 12, 
but plotted here in the plane of the chemical potential and 
the chemical potential difference. Note that the chemical po- 
tential difference is in units of the binding energy, so that 
the diagram is particularly useful for the case with a fixed 
interaction strength, but varying densities. 



5^c,p=i can be determined from the weak coupling ex- 
pression l|5.23p . 



2m 



(5.26) 



as the validity of the equation is restricted to 7 <C 
max{p, 1 — p\. The determination of S^ic,p=o as 7 ^ 
turns out to be very difRcult. Fortunately, it has been 
studied by Krivno y an d Ovchinnikov [sll, and Fuchs, Re- 
cati, and Zwerger [l06l | in detail. Here we only quote their 
result, 
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This predicts the same exponent —Tr'^/{2j) as the BCS 
mean-field theory. However, there is a different power- 
law dependence of the prefactor on the dimensionless cou- 
pling constant, i.e., it has an extra ^ factor. In Fig. 
14, we plot these analytic predictions using dashed and 
dot-dashed lines. They are in excellent agreement with 
the exact numerical results in the regions where they are 
valid. 

For a later reference, in Fig. 15 we reconstruct the 
phase diagram in the plane of the chemical potential and 
the chemical potential difference. Both of them are mea- 
sured in units of the binding energy. It is clear that in the 
strong coupling limit, the two critical chemical potential 
differences converge to the half of the binding energy, and 
the phase space for the FFLO states therefore becomes 
much narrower. 



VI. SELF-CONSISTENT BDG APPROACH IN A 
HARMONIC TRAP 

To make a quantitative contact with the on-going ex- 
periments, it is crucial to take into account the trapping 
potential that is necessary to prevent the atoms from es- 
caping. In this section we turn to describe a ID polarized 
gas in harmonic traps, using the mean-field BdG equa- 
tions. 

With the trap Vtrap {x) = muj'^x^/2, the BdG for- 
malism is essentially the same as that under a peri- 
odic boundary condition, except a few modifications: (1) 
First, one has to replace everywhere the chemical poten- 
tial ^ by a local potential /i — Vtrap{x). (2) Accordingly, 
to solve the BdG equation, it is convenient to use the 
eigenfunctions of the harmonic trap. 



<^„ (x) = AnHn 
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aiio 



exp 



Ho 



(6.1) 



as the set of the expanding basis. Here Hn (x) is the 

1/2 

Hermite polynomial with an order n, a^o — [fi/ {m.uj)] 
the characteristic harmonic oscillator length, and An — 
\/l/(7ri/^2"n!) the normalization factor for single parti- 
cle eigenfunctions. (3) Thirdly, for the convenience of the 
numerical calculations, it is better to take the trap units, 
i.e., m ^ h ~ Lu ^ I, so that the length and energy will 
be measured in units of the characteristic harmonic oscil- 
lator length Uho and hw, respectively. (4) Finally, in the 
presence of the trap, there is no restriction for the initial 
guess of the order parameter. We may then initialize the 
order parameter by choosing some random values. 

We have performed a calculation for a gas with N = 
128 fermions in traps at zero temperature. The Fermi en- 
ergy under the unpolarized condition is Ep — [N/2)huj — 
GAhw. We therefore take a cut-off energy Ec = GEp — 
SSAhuj and keep up to 6A^ = 768 single particle eigen- 
functions. These parameters are already very large to 
ensure the accuracy of the calculations. As mentioned 
earlier, we use the dimensionless coupling parameter at 
the trap center, 70 = Traho/ [N^^'^ciid)-, to characterize 
the interaction. In Fig. 16, we present the BdG results 
for the density profiles (soHd fines) and the order parame- 
ter (dot-dashed lines) at a moderate interaction strength 
7o = 1.6 for three total spin polarizations as indicated. 

For a pure BCS superfiuid with zero polarization (Fig. 
16a), the spin up and down density profiles coincide, and 
decrease monotonically as expected. However, the order 
parameter is non-monotonic: it increases slowly up to 
the boundary of the trap, and then drops to zero very 
rapidly. A maximum at the trap edge then arises in the 
order parameter, in marked contrast to the 3D cases, 
where the order parameter decreases monotonically. This 
maximum is due to the low dimensionality of the gas. 
Recall that the BCS prediction of the gap for a uniform 
gas Abcs — 8ei? exp[— 7r^/(27)]. At the local position 
X, ep oc (x), while 7 = 2/ [aiDn{x)]. As a result, the 
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Figure 16: (Color online) Density profiles (solid lines) and 
order parameters (dot-dashed lines) of a trapped Fermi gas 
at several total spin polarizations as indicated. The dimen- 
sionless coupling constant at the trap center 70 is 1.6. With 
increasing the total spin polarization, the FFLO enters grad- 
ually at center, leading to two phase separation phases. 



position dependent order parameter is given by, 

^2 



^BCS 



(x) oc n (x) exp 



- — air,n(x) 



(6.2) 



which is a product of n'^ (x) and of an exponent. These 
two parts decrease and increase respectively towards the 
trap edge. Particularly, the increase of the exponent is 
due to the increase of the effective interactions, which 
becomes much larger with decreasing density. Therefore 
their interplay should result in a maximum. In general, 
the exponent is dominant, thereby the sharp decrease or 
the maximum of Abcs (x) occurs at the trap edge for a 
moderate local density. 

With increasing total spin polarization, the order pa- 
rameter starts to oscillate at the trap center, suggesting 
the entry of FFLO-type states at center. Correspond- 
ingly, the spin up and down density profiles are no longer 
the same. For a small total spin polarization (Fig. 16b), 
the oscillation of the order parameter is restricted at the 
trap center, and the ordinary BCS order parameter still 
persists at the edge. As a consequence, we find a phase 
separation phase consisting of a FFLO state at the trap 
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Figure 17: (Color online) Local fermionic density of states 
of a trapped polarized Fermi gas at an interaction strength 
70 = 1.6. The remarkable two-energy-gap structure is robust 
in the trap environment. 



center and a standard BCS state outside. There is also 
a very small region with a weak oscillation of the order 
parameter, occurring exactly at the trap boundary. Pre- 
sumably, it is a finite size effect. As we shall see later, 
the resulting normal cloud at the boundary is an artifact 
of the mean-field theory, which turns to break down at 
sufficient small densities or large interactions. 

Increasing further the spin polarizations (Fig. 16c), the 
oscillations of the order parameter penetrate the whole 
cloud. We find then another phase separation phase, with 
an interior core of a FFLO superfiuid phase and an outer 
shell of the normal component. Therefore, there should 
be a critical total spin polarization, Pc, that separates 
the two phase separation phases. The periodicity of the 
oscillations in the FFLO phase can be estimated, and we 
find a reasonable agreement with the single-plane-wave 
estimation for q if we treat the gas as locally homogeneous 
at the trap center. 

The validity of the mean-field BdG calculations in the 
trap environment will be commented later on, by com- 
paring the mean-field density profiles with that obtained 
from the exact Gaudin solution and the local density ap- 
proximation. The physical reason for the two phase sep- 
aration phases and the value of Pc, as well as the small 
oscillations in the density profiles, will also be addressed. 

Finally, we study the local fermionic density of the 
state in the trap. In Fig. 17, we report the density of 
states at the trap center for a BCS superfiuid (a) and a 
FFLO superfiuid (b). In the presence of the trap, we find 
that the essential feature of a two-energy-gap structure 
in the FFLO state is still apparent. This may provide us 
a useful experimental signature to detect indirectly the 
FFLO states. 
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VII. ASYMPTOTICALLY EXACT GAUDIN 
SOLUTIONS IN A HARMONIC TRAP 



For a large number of fermions, a useful method to 
account for an external trapping potential traps is to 
use the local density approximation Together 
with the Gaudin solution for the homogeneous equation 
of states of a polarized Fermi gas, this gives an asymp- 
totically exact result as long as N :s> 1. This condition 
is readily satisfied in the on-going ID experiment, where 
the typical number of atoms N ^ 100. 

The main idea of the local density approximation is 
that the system can be treated locally as infinite matter 
with a local chemical potential. We then partition the 
cloud into many cells in which the number of fermions is 
much greater than unity. Provided that the variation of 
the trap potential across the cell is small compared with 
the local Fermi energy, the interface effects are negligible 
[tbI . [sol . Qualitatively, the interface energy should scale 
like N~^^'^ compared to the total energy, where d is the 
dimensionality. 

In detail, the local density approximation amounts to 
determining the chemical potential fig = [jii^g -\- f^ig)/2 
and the chemical potential difference S^g = {ji^g— ji^g) /2 
of the inhomogeneous gas from the local equilibrium con- 
ditions, 



and the normalization conditions, 
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where n{x) and p{x) are respectively the total linear den- 
sity and the local spin polarization, and P the total spin 
polarization. We have used a subscript "g" to denote the 
global chemical potentials. 

To solve these equations, we rewrite the chemical po- 
tentials in the form, 



^^^[n{x),p{x)] = —n^{x)^i^[-i{x),p{x)], (7.5) 
tJ'i[n{x),p{x)] = —n^{x)^iiY^{x),p{x)], (7.6) 



where Ji^ are the reduced chemical potentials, depending 
on the dimensionless coupHng constant and local spin 
polarization only. Further, it is convenient to rescale the 
chemical potentials, coordinate and total linear density 



into a dimensionless form, i.e., 
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Then the local equilibrium equations and the normaliza- 
tion equations can be rewritten as. 
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where 7(2;) — 2/h{x). The terms on the left-side hand 
of the last two equations emphasize that the properties 
of the cloud rely on two dimensionless parameters, 70 
and P. In particular, the coupling constant in a trap 
is controlled by 70, where 70 ^ 1 corresponds to weak 
coupHng, while 70 ^ 1 corresponds to the strongly inter- 
acting regime. 

The numerical procedure for the local density approx- 
imation is straightforward. For given parameters 70 and 
P, and initial guess for flag, we invert the dimensionless 
local equilibrium equations to find 7(0;) and p{x). The 
chemical potentials fia-g are then adjusted slightly to en- 
force number conservation, giving a better estimate for 
the next iterative step. The iteration is continued un- 
til the number conditions are satisfied within a certain 
range. 



A. Density profiles: LDA vs BdG 

In Fig. 18, we give the density profiles obtained from 
the local density approximation using dashed lines. For 
comparison, we show also the results of the BdG so- 
lutions. Apart from a negligible difference at the trap 
boundary (due to a breakdown of mean-field theory) , we 
find a good agreement. This becomes even better as the 
total spin polarization increases. In particular, the two 
phase separation phases found in the BdG calculations 
are evident. 

The appearance of the phase separation phases is easy 
to understand. Within the local density approximation, 
the local chemical potential fi{x) decreases parabolically 
away from the center of the trap while the local chemical 
potential difference 5fi{x) stays constant. It is then evi- 
dent from Fig. 15 that with a nonzero spin polarization 
we always have a polarized FFLO superfiuid at the trap 
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Figure 18: (Color online) Density profiles of a trapped gas, 
calculated by the exact Gaudin solution and the local den- 
sity approximation, are shown for several spin polarizations 
as indicated. For comparison, we plot also the self-consistent 
mean-field BdG predictions. They are in reasonable agree- 
ment at the center. A discrepancy occurs at the trap edge, 
where for small polarization, the approximate BdG calcula- 
tion overestimates the size of the unpolarized BCS shell. 

center where the local chemical potential (or interaction 
parameter) is large (or small). Away from the center 
with decreasing local chemical potential, the gas enters 
into either an unpolarized BCS superfluid or a fully po- 
larized normal cloud, depending whether the chemical 
potential difference is smaller than a half of the binding 
energy or not. Thus, there is a critical chemical potential 
difference S^j,c = eh/2 that separates the inhomogeneous 
system into two phase separation states: a mixture of a 
polarized superfluid core and an unpolarized superfluid 
shell (FFLO-BCS), or a coexistence of a polarized su- 
perfluid at the center and a fully polarized normal gas 
outside (FFLO-BCS). 

It should be noted that the former phase separation 
phase is exotic, as the BCS-Hke superfluid state occurs 
at the edge of the trap, in marked contrast to the 3D 
case. This is caused by the peculiar effects of low dimen- 
sionality, for which the gas becomes more nonideal with 
decreasing ID density towards the edge of the trap, and 
hence the energy required to break the pairs approaches 
£6/2 from below. As S^g < eb/2, there should be a fully 




Figure 19: (Color online) Phase diagram of a one-dimensional 
trapped spin-polarized Fermi gas. The dashed line and dot- 
dashed line are the asymptotic results for the critical spin 
polarization in the strongly and weakly interacting regimes, 
respectively. 



paired region once the local critical chemical potential 
Sf^c,p=o > Sfig, i.e., the BCS-Hke superfluid. 

Though the basic feature of the BdG results is well 
reproduced by the local density approximation calcula- 
tions, we note that there are still some discrepancies that 
merit careful examination. First, with decreasing the 
density the mean-fleld theory seems to fail at the trap 
edge, as shown in Figs. (18a) and (18b). For a small 
polarization P = 0.05 (Fig. 18b), a notable discrepancy 
thus occurs at the trap edge. The very small unpolarized 
BCS shell, roughly from 0.80N^/^aho to 0.85N^/^aho as 
predicted by the LDA calculation, becomes strongly over- 
estimated by the mean-fleld calculation. Secondly, there 
are small oscillations in the BdG density proflles. Pre- 
sumably, these oscillations, observed also in a box with 
periodic boundary conditions, are either due to the pres- 
ence of the FFLO states or due to a flnite size effect. 
Considering the absence of the true long-range order in 
ID, we prefer the later interpretation, and regard them 
as the Priedel oscillations caused by the residual unpaired 
atoms. To check this point, in the BdG calculations we 
have varied the total number of fermions, while keep- 
ing other parameters invariant. The oscillations become 
less pronounced with increasing numbers of atoms. We 
emphasize that in the on-going experiments, the total 
number of atoms is about one hundred. Therefore, the 
oscillations in the density proflles could be observed ex- 
perimentally. However, they may hardly be considered 
as a fundamental signature of the presence of the FFLO 
states. 



B. Phase diagram of a polarized Fermi gas in traps 

We may determine numerically the critical spin polar- 
ization Pc from the critical chemical potential difference 
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S/ic = €6/2. In Fig. 19, we present Pc as a function of the 
interaction coupling constant 70, giving rise to a phase 
diagram of the inhomogeneous polarized ID Fermi gas 
[ssl . [9^ . Again, the asymptotic behavior of Pc may be 
computed analytically in the weak and strong coupling 
limits. These are shown in the figure using a dashed Hne 
and a dot-dashed Hne, respectively. 

Consider first a strongly interacting gas with j{x) > 
7o ^ 1. Using the asymptotic expression for the chemical 
potential and chemical potential difference, the rescaled 
local equilibrium equations can be rewritten as. 



TT^n^ (x) , , , ,2 
—-^[l-p{x)f + - ^ fig, 



2 + ^^MPm + ^^B[pix)] - MJ,14) 

where = -I + 2p (x) + 15p^ (x) and B [p (x)] = 

-p (x) /4 + (x) /2 - 67p^ (x) /12. Note that in this 
limit n (x) ^ 1. In the rescaled units, the critical chem- 
ical potential difference 6pg is exactly 1/2. Therefore, if 
we consider up to y4[p(x)] only in the expansion, we find 
that the local spin polarization should satisfy. 



15/ (x) + 2p{x)-l^ 0, 



(7.15) 



which yields p{x) = 1/5 and hence the total spin polar- 
ization Pc = 1/5. The improvement to the next order 
requires the inclusion of the term B [p (x)]. For this pur- 
pose, we assume p{x) = 1/5 — S{x), where S{x) <^ 1. 
The summation of A [p {x)] and B [p (x)] terms should be 
zero at the critical polarization. Thus, to leading order 
of S{x), we find that. 



Six) 
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n (x) 



(7.16) 



The density profile n {x) can be determined by using the 
local equilibrium equation for fig, which to a good ap- 
proximation 



1 TT^ o , , x'^ 

- H (x) -\ 

2 50 ^ ^ 2 



Ms- 



Combined with the normalization 



n{x)dx = l/(7r^7Q), we find that. 



n (x) 



71"^ 70 L 



1 



bn^jQX^ 



-,1/2 



(7.17) 
condition, 

(7.18) 



Therefore, we determine the critical spin polarization us- 
ing Pc — 7r^7o n (x) p (x) dx and find that. 
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(7.19) 



The consideration in the weak coupling limit is much 
simple. In the rescaled units. 



6n [n (x) ,pix)] 



TT n (x) p (x) 



SUg 



(7.20) 



where in this limit n {x) ^ 1. By setting 5^g = 1/2, we 
then obtain. 



p{x) 



2 1 

TT^ (x) 



(7.21) 



Using again the normalization condition for the total 
number of atoms, the rescaled (ideal) density profile takes 
the form, 



- /'-^ ^ fi 2 2-21 1/2 



(7.22) 



Thus, by integrating out Pc = tt^Jq 2/['K'^n{x)]dx, 
we find that, 



p _ 7o 



(7.23) 



VIII. CONCLUSIONS AND SOME REMARKS 

In conclusion, we have presented a systematic study of 
an attractive polarized atomic Fermi gas in one dimen- 
sion, both in free space and in a harmonic trap. The 
theoretical approaches include the (asymptotically) ex- 
act Bethe ansatz solution and two mean-field approx- 
imations: the single-plane- wave approximation for the 
order parameter and the self-consistent Bogoliubov-de 
Gennes equations. These useful tools provide us with 
quantitative phase diagrams in both uniform and har- 
monic trapped systems. Our main results may be sum- 
marized as follows, in response to the theoretical issues 
raised in the Introduction: 

(A) We have clarified the structure of the one- 
dimensional FFLO states in a uniform gas. For small 
spin polarization, the FFLO order parameter behaves like 
a lattice of instantons and anti-instantons, which carry 
the excess unpaired atoms. For a large spin polarization, 
the singularity of the instantons merges together. Thus, 
the form of the order parameter becomes a cosine func- 
tion, as originally proposed by Larkin and Ovchinnikov 
[29| . The nodes in the FFLO order parameter lead to a 
two-energy-gap structure in the local fermionic density 
of states, which may be experimentally observable using 
spectroscopic methods. 

(B) We have determined the nature of the phase tran- 
sition from a BCS superfiuid state to a FFLO phase. It 
is a smooth second order transition. As a consequence, 
a one-dimensional phase separation does not occur for 
a homogeneous gas. Turning to the trapped case, we 
find two exotic phase separation phases. However, these 
phase separations are simply trap effects. 

(C) We have checked the validity of the two mean- 
field approaches in the weakly or moderately interact- 
ing regimes, by comparing the results with the exact or 
asymptotically exact Bethe ansatz solutions. The mean- 
field methods are found to provide a useful description 
in these regimes. In particular, by comparing the equa- 
tions of state and density profiles, we have shown that 
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the spin polarized superfluid in the Bethe ansatz solution 
corresponds to an FFLO state, with a real (cosine-like) 
order parameter. This correspondence, however, does 
not hold quantitatively in the strongly interacting regime. 
The Bethe ansatz solutions do not result in any abrupt 
changes for the polarized superfluid, as the interaction 
strengths increase from the weak to strong regimes. 

Though our study is restricted here to the one- 
dimensional case, we can still obtain some insight into 
the phase diagram of a three-dimensional polarized Fermi 
gas. This is under strong debate at the moment. Two 
remarks may be in order in this respect. 

One key remark is that the FFLO window in three di- 
mension can be expected to be much larger than that ob- 
tained from mean-field calculations with a single-plane- 
wave assumption for the order parameter. As we have 
noted, by improving the form of the order parameter 
to the Larkin and Ovchinnikov (LO) type, A(x) cx 
cos[q-x], Yoshida and Yip have indeed found recently 
that the FFLO state becomes more stable [H^l- Further, 
the one-dimensional results indicate that one might ex- 
pect a smooth phase transition from the BCS state to 
FFLO state in three dimensions, although clearly this 
needs to be checked with a full three-dimensional calcu- 
lation. 

Another interesting issue concerns the existence of a 
phase separation in a three dimensional homogeneous po- 
larized gas. From the one-dimensional calculations, we do 



not find any strong indication for this. Accordingly, the 
experimentally observed phase separation may simply be 
understood as a trap effect. We note, however, that the 
three dimensional strongly interacting BEC limit has no 
correspo ndence in the one-dimensional attractive polar- 
ized gas |l05l . Il06| . It that limit, a homogeneous polar- 
ized superfluid, which may be called the Sarma phase, be- 
comes stable [33, 34, 44]. This phase has a different sym- 
metry from the spatially inhomogeneous FFLO phase. 
Therefore, there could be another phase intervening be- 
tween the Sarma phase and the FFLO phase. This may 
be a possible reason for the observation of phase separa- 
tion in three dimension. If this exists, we would expect 
that phase separation for a homogeneous gas would be 
restricted to the strongly-interacting regime near unitar- 
ity. 
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